pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.     

Simple method for the quantitative examination of extra column band broadening in microchromatographic systems

In recent years capillary chromatography has gained popularity for trace analyses. Most often UV or electrochemical detection is employed because the small peak volumes make post-column derivatization challenging. We have developed a simple method based on flow injection for determining contributions to peak broadening from post-column reactors. The only requirement for application of our methodology is that diffusion be in the Taylor regime so that radial concentration gradients are relaxed enabling mixing purely by diffusion.     

On the viability of small endohedral hydrocarbon cage complexes: X@C4H4, X@C8H8, X@C8H14, X@C10H16, X@C12H12, AND X@C16H16

Small hydrocarbon complexes (X@cage) incorporating cage-centered endohedral atoms and ions (X = H(+), H, He, Ne, Ar, Li(0,+), Be(0,+,2+), Na(0,+), Mg(0,+,2+)) have been studied at the B3LYP/6-31G(d) hybrid HF/DFT level of theory. No tetrahedrane (C(4)H(4), T(d)()) endohedral complexes are minima, not even with the very small hydrogen atom or beryllium dication. Cubane (C(8)H(8), O(h)()) and bicyclo[2.2.2]octane (C(8)H(14), D(3)(h)()) minima are limited to encapsulating species smaller than Ne and Na(+). Despite its intermediate size, adamantane (C(10)H(16), T(d)()) can enclose a wide variety of endohedral atoms and ions including H, He, Ne, Li(0,+), Be(0,+,2+), Na(0,+), and Mg(2+). In contrast, the truncated tetrahedrane (C(12)H(12), T(d)()) encapsulates fewer species, while the D(4)(d)() symmetric C(16)H(16) hydrocarbon cage (see Table of Contents graphic) encapsulates all but the larger Be, Mg, and Mg(+) species. The host cages have more compact geometries when metal atoms, rather than cations, are inside. This is due to electron donation from the endohedral metals into C-C bonding and C-H antibonding cage molecular orbitals. The relative stabilities of endohedral minima are evaluated by comparing their energies (E(endo)) to the sum of their isolated components (E(inc) = E(endo) - E(cage) - E(x)) and to their exohedral isomer energies (E(isom) = E(endo) - E(exo)). Although exohedral binding is preferred to endohedral encapsulation without exception (i.e., E(isom) is always exothermic), Be(2+)@C(10)H(16) (T(d)(); -235.5 kcal/mol), Li(+)@C(12)H(12) (T(d)(); 50.2 kcal/mol), Be(2+)@C(12)H(12) (T(d)(); -181.2 kcal/mol), Mg(2+)@C(12)H(12) (T(d)(); -45.0 kcal/mol), Li(+)@C(16)H(16) (D(4)(d)(); 13.3 kcal/mol), Be(+)@C(16)H(16) (C(4)(v)(); 31.8 kcal/mol), Be(2+)@C(16)H(16) (D(4)(d)(); -239.2 kcal/mol), and Mg(2+)@C(16)H(16) (D(4)(d)(); -37.7 kcal/mol) are relatively stable as compared to experimentally known He@C(20)H(20) (I(h)()), which has an E(inc) = 37.9 kcal/mol and E(isom) = -35.4 kcal/mol. Overall, endohedral cage complexes with low parent cage strain energies, large cage internal cavity volumes, and a small, highly charged guest species are the most viable synthetic targets.     

Homonuclear transition-metal trimers

Density-functional theory has been used to determine the ground-state geometries and electronic states for homonuclear transition-metal trimers constrained to equilateral triangle geometries. This represents the first application of consistent theoretical methods to all of the ten 3d block transition-metal trimers, from scandium to zinc. A search of the potential surfaces yields the following electronic ground states and bond lengths: Sc3(2A1',2.83 A), Ti3(7E',2.32 A), V3(2E",2.06 A), Cr3(17E',2.92 A), Mn3(16A2',2.73 A), Fe3(11E",2.24 A), Co3(6E",2.18 A), Ni3(3A2",2.23 A), Cu3(2E',2.37 A), and Zn3(1A1',2.93 A). Vibrational frequencies, several low-lying electronic states, and trends in bond lengths and atomization energies are discussed. The predicted dissociation energies DeltaE(M3-->M2+M) are 49.4 kcal mol(-1)(Sc3), 64.3 kcal mol(-1)(Ti3), 60.7 kcal mol(-1)(V3), 11.5 kcal mol(-1)(Cr3), 32.4 kcal mol(-1)(Mn3), 61.5 kcal mol(-1)(Fe3), 78.0 kcal mol(-1)(Co3), 86.1 kcal mol(-1)(Ni3), 26.8 kcal mol(-1)(Cu3), and 4.5 kcal mol(-1)(Zn3).     

Computational study of the aminolysis of esters. The reaction of methylformate with ammonia

The aminolysis of esters is a basic organic reaction considered as a model for the interaction of carbonyl group with nucleophiles. In the present computational study the different possible mechanistic pathways of the reaction are reinvestigated by applying higher level electronic structure theory, examining the general base catalysis by the nucleophile, and a more comprehensive study the solvent effect. Both the ab initio QCISD/6-31(d,p) method and density functional theory at the B3LYP/6-31G(d) level were employed to calculate the reaction pathways for the simplest model aminolysis reaction between methylformate and ammonia. Solvent effects were assessed by the PCM method. The results show that in the case of noncatalyzed aminolysis the addition/elimination stepwise mechanism involving two transition states and the concerted mechanism have very similar activation energies. However, in the case of catalyzed aminolysis by a second ammonia molecule the stepwise mechanism has a distinctly lower activation energy. All transition states in the catalyzed aminolysis are 10-17 kcal/mol lower than those for the uncatalyzed process.     

Cyclic perfluorocarbon radicals and anions having high global warming potentials (GWPs): structures, electron affinities, and vibrational frequencies

Adiabatic electron affinities, optimized molecular geometries, and IR-active vibrational frequencies have been predicted for small cyclic hydrocarbon radicals C(n)H(2)(n)(-)(1) (n = 3-6) and their perfluoro counterparts C(n)F(2)(n)(-)(1) (n = 3-6). Total energies and optimized geometries of the radicals and corresponding anions have been obtained using carefully calibrated (Chem. Rev. 2002, 102, 231) density functional methods, namely, the B3LYP, BLYP, and BP86 functionals in conjunction with the DZP++ basis set. The predicted electron affinities show that only the cyclopropyl radical tends to bind electrons among the hydrocarbon radicals studied. The trend for the perfluorocarbon (PFC) radicals is quite different. The electron affinities increase with expanding ring size until n = 5 and then slightly decrease at n = 6. Predicted electron affinities of the hydrocarbon radicals using the B3LYP hybrid functional are 0.24 eV (C(3)H(5)/C(3)H(5)(-)), -0.19 eV (C(4)H(7)/C(4)H(7)(-)), -0.15 eV (C(5)H(9)/C(5)H(9)(-)), and -0.11 eV (C(6)H(11)/C(6)H(11)(-)). Analogous electron affinities of the perflurocarbon radicals are 2.81 eV (C(3)F(5)/C(3)F(5)(-)), 3.18 eV (C(4)F(7)/C(4)F(7)(-)), 3.34 eV (C(5)F(9)/C(5)F(9)(-)), and 3.21 eV (C(6)F(11)/C(6)F(11)(-)).     

Hypervalency avoided: simple substituted BrF3 and BrF5 molecules. Structures, thermochemistry, and electron affinities of the bromine hydrogen fluorides HBrF2 and HBrF4

Five different pure density functional theory (DFT) and hybrid Hartree-Fock/DFT methods have been used to search for the molecular structures, thermochemistry, and electron affinities of the bromine hydrogen fluorides HBrF(n)/HBrF(n)(-) (n = 2, 4). The basis sets used in this work are of double-zeta plus polarization quality in conjunction with s- and p-type diffuse functions, labeled as DZP++. Structures with Br-F and Br-H normal bonds, that is, HBrF(2)/HBrF(2)(-) with C(2v) or C(s) symmetry and HBrF(4)/HBrF(4)(-) with C(4v) or C(s) symmetry, are genuine minima. However, unlike the original BrF(3) and BrF(5) molecules, the global minima for HBrF(n)/HBrF(n)(-) (n = 2, 4) species are predicted to be complexes, some of which contain hydrogen bonds. The demise of the hypervalent structures is due to the availability of favorable dissociation products involving HF, which has a much larger dissociation energy than F(2). Similar reasoning suggests that PF(4)H, SF(3)H, SF(5)H, ClF(2)H, ClF(4)H, AsF(4)H, SeF(3)H, and SeF(5)H will all be hydrogen bond structures incorporating diatomic HF. The most reasonable theoretical values of the adiabatic electron affinities (EA(ad)) are 3.69 (HBrF(2)) and 4.38 eV (HBrF(4)) with the BHLYP method. These electron affinities are comparable to those of the analogous molecules: Br(2)F(n), ClBrF(n), and BrF(n)(+1) systems. The first F-atom dissociation energies for the neutral global minima are 60 (HBrF(2)) and 49 kcal/mol (HBrF(4)) with the B3LYP method. The first H-atom dissociation energies for the same systems are 109 (HBrF(2)) and 116 kcal/mol (HBrF(4)). The large Br-H bond energies are not sufficient to render the hypervalent structures energetically tenable. The dissociation energies for the complexes to their fragments are relatively small.     

The solid state structure of arene-mercury(II) complexes from magic-angle spinning carbon-13 NMR and the solution carbon-13 NMR of some complexes of mercury(II) with arenes having bulky aliphatic substituents

The cross-polarization magic angle spinning ¹³C NMR spectra of Hg(SbF₆)₂ - 2 Arene (Arene = C₆HMe₅, 1,2,4,5-C₆H₂Me₄, 1,2,3,4-C₆H₂Me₄, or C₆H₆) have been measured. The spectra of the complexes of C₆HMe₅ and 1,2,4,5-C₆H₂Me₄ are consistent with static η¹-bonding of the mercury to the arene at an unsubstituted carbon atom, while the spectra of the 1,2,3,4-C₆H₂Me₄ and C₆H₆ complexes show the arene to have time-averaged C_s or C₂, and C₆ symmetry respectively, at the temperature of measurement (300 K).The reduced temperature ¹³C NMR spectra of Hg(Arene)_n²⁺ (n = 1 or 2; Arene = 1,3,5-C₆H₃R₃ (R = Me, i-Pr, or t-Bu)) in SO₂ solution are also reported and affirm that in these intramolecularly mobile species the mercury bonds in an η¹-manner, with unsubstituted aryl carbon atoms being the strongly preferred point of mercury attachment. This site preference is further demonstrated by the solution ¹³C NMR spectra of Hg(Arene)_n²⁺ (Arene = 1,2,3,4-C₆H₂-Me₄, n = 1 or 2; Arene = 1,4-C₆H₄R₂, R = Me or t-Bu, n = 1). The spectra of the 1,4-C₆H₄R₂ complexes and Hg(p-C₆H₄-t-BuMe)²⁺ provide clear evidence for steric influence of the binding site.Like Hg(C₆Me₆)₂²⁺, but unlike most of the complexes of substituted benzenes which have been studied, Hg(1,3,5-C₆H₃-i-Pr₃)₂²⁺ exchanges only slowly with excess free ligand.     

What is the nature of polyacetylene neutral and anionic chains HC(2n)H and HC(2n)H(-) (n = 6-12) that have recently been observed?

The optimized geometries, adiabatic electron affinities, and IR-active vibrational frequencies have been predicted for the long linear carbon chains HC(2n)H. The B3LYP density functional combined with the DZP basis set was used in this theoretical study. The computed physical properties are discussed. The predicted electron affinities form a remarkably regular sequence: 1.78 (HC(12)H), 2.08 (HC(14)H), 2.32 (HC(16)H), 2.53 (HC(18)H), 2.69 (HC(20)H), 2.83 (HC(22)H), and 2.95 eV (HC(24)H). The predicted structures display an alternating triple and very short single bond pattern, with the degree of bond alternation significantly less for the radical anions.